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ABSTRACT 

The Earth appears non-chondritic in its abundances of refractory lithophile elements, posing a 
significant problem for our understanding of its formation and evolution. It has been suggested that 
this non-chondritic composition may be explained by collisional erosion of differentiated planetesimals 
of originally chondritic composition. In this work, we present V-body simulations of terrestrial planet 
formation that track the growth of planetary embryos from planetesimals. We simulate evolution 
through the runaway and oligarchic growth phases under the Grand Tack model and in the absence 
of giant planets. These simulations include a state-of-the-art collision model which allows multiple 
collision outcomes, such as accretion, erosion, and bouncing events, that enables tracking of the 
evolving core mass fraction of accreting planetesimals. We show that the embryos grown during this 
intermediate stage of planet formation exhibit a range of core mass fractions, and that with significant 
dynamical excitation, enough mantle can be stripped from growing embryos to account for the Earth’s 
non-chondritic Ee/Mg ratio. We also find that there is a large diversity in the composition of remnant 
planetesimals, with both iron-rich and silicate-rich fragments produced via collisions. 

Subject headings: Earth — methods: numerical — planets and satellites: composition — planets and 
satellites: formation — planets and satellites: terrestrial planets — protoplanetary 
disks 


1. INTRODUCTION 

As the best studied planetary system, our own solar 
system provides the most important benchmark for theo¬ 
ries of planet formation. Critically, we have detailed com¬ 
positional information on the Solar system from spectro¬ 
scopic measurements of the sun and analyses of mete¬ 
orites (e.g. Lodders 2003). By assuming that planets 
have the composition of primitive meteorites, at least for 
refractory, lithophile elements that vary little between 
meteorite groups, inferences about planetary evolution 
can be made from the deviation of observed compositions 
from such a likely starting point. In the last decade, the 
validity of the notion of a chondritic bulk planetary com¬ 
position, has faced close scrutiny (e.g. Boyet & Carlson 
2005) and several studies have suggested that the pro¬ 
cess of accretion may inevitably generate planets with 
non-chondritic compositions (e.g. Bourdon et al. 2008; 
O’Neill & Palme 2008). However, these ideas have yet to 
be explored in detail with physical models. 

The terrestrial planets of our solar system exhibit a 
range of Ee/Mg ratios. The Earth (Ee/Mg 2.1, O’Neill 
& Palme 2008) is significantly enriched in iron compared 
to both carbonaceous chondrites (1.92±0.08, Palme & 
Jones 2003) and the solar photosphere (1.87±0.4, Lod¬ 
ders 2003), and Mercury with its massive core must have 
a large iron excess (estimated Ee/Mg = 9.92, Morgan & 
Anders 1980). Estimates for Venus (Morgan & Anders 
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1980) and Mars (Lodders & Eegley 1997) suggest they 
have Ee/Mg ratios similar to the Earth and chondrites 
respectively. In contrast the Moon is expected to have 
a Ee/Mg ratio significantly lower than chondritic values, 
and Ee/Mg ratios for meteorites (including differentiated 
bodies) range from less than 0.08 to greater than 25 (Nit- 
tler et al. 2004). Such extremes would not be expected for 
primitive bodies; iron meteorites, for example, probably 
represent the surviving cores of differentiated planetesi¬ 
mals. It has been proposed that stripping of crust and 
mantle via collisions during accretion led to this observed 
compositional variation (e.g. Benz et al. 1988; O’Neill & 
Palme 2008). 

Recently some progress has been made by use of N- 
body accretion models which allow for imperfect merging 
(Bonsor et al. 2015; Dwyer et al. 2015). Both of these 
studies investigated discs unperturbed by giant planets, 
the calmest possible scenario for the formation of terres¬ 
trial planets. Migration of the giant planets has been 
argued to play a key role in explaining the small size of 
Mars, depopulating the asteroid belt and scattering po¬ 
tentially volatile rich material into the inner Solar system 
(e.g. Walsh et al. 2011). The excitation of planetesimals 
such migration causes may also have significant implica¬ 
tions for the importance of collisional erosion during the 
accretion of the terrestrial planets. 

In this work, we extend the work of Bonsor et al. (2015) 
(hereafter Paper I) in assessing the likelihood of non- 
chondritic accretion of terrestrial planets, by examining 
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the dynaniically hot scenario provided by the Grand Tack 
model (Walsh et al. 2011). We compare this to an up¬ 
dated model of calm disc accretion, which additionally 
explores the effects of gas drag. 

2. NUMERICAL METHOD 

We use a modified version of the parallelized Wbody 
code PKDGRAV (Richardson et al. 2000; Stadel 2001) to 
simulate the gravitational interactions of planetesimals in 
the potential well of a star. A state-of-the-art collision 
model (EDAGM, Leinhardt & Stewart 2012; Leinhardt 
et al. 2015) has been incorporated into this code in or¬ 
der to calculate the results of various collision types, in¬ 
cluding partial accretion, hit and run and erosive events. 
PKDGRAV allows us to simulate much larger numbers 
of particles than codes used in similar works (e.g. Walsh 
et al. 2011; Ghambers 2013; O’Brien et al. 2014), and 
to include planetesimal-planetesimal interactions. How¬ 
ever, the second-order integrator would not achieve suf¬ 
ficient accuracy for the long term orbital evolution of 
embryos during the giant impact phase - the final stage 
of planet formation characterised by stochastic collisions 
between large embryos that have developed chaotic orbits 
due to their mutual interaction. In this work we simulate 
the intermediate stages of planet formation characterised 
by runaway and oligarchic growth of planetesimals and 
embryos. All simulations were carried out using the Uni¬ 
versity of Bristol’s BlueGrystal^ supercomputer. 

2.1. Initial planetesimal disc 

In this work we present simulations initiated with ei¬ 
ther 10 000 (low resolution) or 100 000 (high resolution) 
particles arranged in an annulus, extending from 0.5 AU 
to either 1.5 or 3 AU, around a one solar mass star. 

The initial disc of planetesimals is distributed accord¬ 
ing to a surface density, E = similar to that 

of the Minimum Mass Solar Nebula (MMSN, Weiden- 
schilling 1977; Hayashi 1981), using a surface density at 
1 AU of El = 10gcm“^ to match previous work (e.g. 
Kokubo & Ida 2002; Leinhardt et al. 2009; Paper I). 

We begin our simulations with planetesimals of 15 dif¬ 
ferent sizes, evenly spaced in mass, m, and chosen ac¬ 
cording to the cumulative number distribution: 

dN = m-Pdm, (1) 

with p = 2.5, a slope appropriate for the runaway growth 
phase (Kokubo & Ida 1996). This initial mass distribu¬ 
tion (black line in Figure 1) is similar to that used in 
Paper I (in which the size distribution was generated by 
running PKDGRAV using perfect merging with a disc 
with five times as many equal size planetesimals). The 
simulations begin with most of the planetesimals having 
radii of a few hundred km in accordance with the mass 
distribution shown in Figure 1. 

In this work we study the intermediate stage of planet 
formation, sizeable planetesimals have already formed 
and will continue to grow into embryos via runaway 
accretion during the simulations. As it is known that 
planetesimals can differentiate in the first few million 
years of the Solar System’s evolution (e.g. Lugmair & 
Shukolyukov 1998; Srinivasan et al. 1999; Amelin 2008; 

^ http://www.bris.ac.uk/acrc/ 


Kruijer et al. 2014) we assume that all planetesimals are 
differentiated at the start of our simulation, and assign 
the same core fraction to all particles, regardless of their 
location within the disc. This uniform initial differentia¬ 
tion is probably a simplification, but will aid understand¬ 
ing of the role of collisions on the compositional evolution 
of growing bodies. 

The initial conditions used for the simulations pre¬ 
sented in this paper are summarised in Tables 1 and 2. 


2.2. Expansion factor 

V-body simulations of the intermediate phase of planet 
formation with large N have often employed an expan¬ 
sion of the particle radii in order to increase the collision 
rate, and hence speed up the evolution, e.g. Kokubo & 
Ida (1996, 2002); Leinhardt & Richardson (2005); Mor- 
ishima et al. (2008). This expansion factor was discussed 
in detail by Kokubo & Ida (1996, 2002), in which they de¬ 
termine that it has very little effect on the evolution dur¬ 
ing this phase, merely altering the timescale for planet 
formation. In this work we are studying the combined 
consequences of many collisions of different types, includ¬ 
ing erosive collisions, and since the expansion factor pre¬ 
vents very close approaches from occurring (which would 
excite planetesimal velocities), it is expected that our 
simulations slightly underestimate mantle erosion. 

While the number of bodies in the simulation is large, 
dynamical friction from the planetesimals controls the ve¬ 
locity dispersion, and the expansion factor does not alter 
their growth mode (Kokubo & Ida 1996). As the number 
of planetesimals decreases, and massive embryos begin to 
significantly alter velocities via gravitational scattering, 
the expansion factor likely causes the results to become 
less accurate. Once the giant impact stage of planet for¬ 
mation begins this becomes a significant problem, and so 
we cannot integrate accurately through the late stages of 
formation using the expansion factor. Due to limitations 
in computing power it is still not practical to conduct 
large N simulations covering millions of years of planet 
formation without using this expansion factor. 

In previous work we have used an expansion factor, 
/ = 6 (Leinhardt & Richardson 2005; Leinhardt et al. 
2009; Paper I; Leinhardt et al. 2015), and, for consis¬ 
tency, we adopt this same value here. In our simulations 
giant planets do not have their radii inflated (any plan¬ 
etesimals colliding with giant planets are forced to merge 
to avoid complications), and gas drag is calculated using 
the natural radius, not the expanded value. 

During the intermediate stages of planet formation we 
are investigating, the expansion factor causes the evolu¬ 
tion to speed up by a factor ^ /^. We therefore quote 
both a simulation time, and an effective time correspond¬ 
ing to realistic radii (/ = 1). 

2.3. Collision model and debris 

EDAGM (Leinhardt & Stewart 2012; Leinhardt et al. 
2015) provides an analytical description of the outcome 
of a collision between gravity-dominated bodies, like the 
planetesimals modelled in this work. Using this colli¬ 
sion model the size and velocity distributions of the frag¬ 
ments of each collision are calculated based on the im¬ 
pact parameter and collision velocity, with each collision 
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Table 1 

Initial conditions for the calm disc simulations. 


Parameter 

Value 

Notes 

Mass of planetesimal disc 

2.5 M0 


Radial extent of disc 

0.5-I.5AU 


Initial number of particles 

100 000 

High resolution 


10 000 

Low resolution 

Initial planetesimal mass 

3.2 X 10-11 M© - 1.5 X 10-® M© 

High resolution 


3.2 X 10-10 M© - 3.6 X 10-® M© 

Low resolution 

Initial planetesimal radius 

196 km - 1530 km 

High resolution (/ = 1) 


423 km - 2050 km 

Low resolution (/ = 1) 

Initial core fraction 

0.22, 0.35 


Planetesimal density 

2 gcm-o 


Resolution limit 

3 X 10-11 M© (1 X 10-^ M©) 

High resolution 


3 X 10-10 M© (1 X 10-^ M©) 

Low resolution 

Timestep 

0.01 yr, 0.005 yr (with gas drag) 


Duration 

600 000 yr (~20Myr effective time) 


Initial gas density, pg,i 

1.4 X lO-o gcm-o 

MMSN 

Gas dissipation start time 

56 kyr (2 Myr effective time) 


Gas dissipation timescale 

2.8kyr (100kyr effective time), oo 



Table 2 

Initial conditions for the Grand Tack simulations. 


Parameter 

Value 

Notes 

Mass of planetesimal disc 

4.85 M© 


Radial extent of disc 

0.5-3.0AU 


Initial number of particles 

100 000 

High resolution 


10 000 

Low resolution 

Initial planetesimal mass 

6.2 X 10-11 M© - 2.9 X 10-® M© 

High resolution 


6.2 X 10-10 M© - 7 X 10-® M© 

Low resolution 

Initial planetesimal radius 

245 km - 1900 km 

High resolution (/ = 1) 


528 km - 2550 km 

Low resolution (/ = 1) 

Initial core fraction 

0.22, 0.35 


Planetesimal density 

2gcm-o 


Resolution limit 

6 X 10-11 M© (2 X 10-^ M©) 

High resolution 


6 X 10-10 M© (2 X 10-4 M©) 

Low resolution 

Timestep 

0.01 yr, 0.005 yr (with gas drag) 


Duration 

600 000 yr (~20Myr effective time) 


Jupiter initial location 

3.5 AU 


Jupiter minimum semi-major axis 

1.5 AU 


Jupiter final location 

5.2 AU 


Migration timescale 

2.8kyr (100kyr effective time), 17kyr (600kyr) 

Fast, slow migration 

Migration start time 

56kyr (2Myr effective time), 220kyr (8Myr) 

Early, late migration 

Initial gas density, pg 

1.4 X lO-o gcm-o 

MMSN 


see equations (6) & (7) 

Based on hydrodynamic simulations 



Mass (M^) 


Figure 1. Evolution of the planetesimal mass distribution (cumulative number of bodies) for a high resolution simulation under the Grand 
Tack scenario (simulation 21, left) and a calm disc simulation (5, right). 
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classified as having one of seven different outcomes: per¬ 
fect merging, partial accretion, hit and run (with the 
projectile intact, disrupted, or supercatastrophically dis¬ 
rupted), erosion, or supercatastrophic disruption. 

Since EDACM can produce many fragments from the 
original two impactors it is necessary to impose a resolu¬ 
tion limit, a mass below which no further resolved par¬ 
ticles will be generated. Any remnants smaller than this 
limit are treated semi-analytically as ‘unresolved debris’. 
This debris is distributed in a series of circular annuli 
each 0.1 AU wide, with the debris placed in the annulus 
corresponding to the location of the collision that pro¬ 
duced it. This material is assumed to have circular Ke- 
plerian orbits. As well as being generated by collisions, 
this unresolved debris is reaccreted by the resolved parti¬ 
cles according to their geometric cross-section (calculated 
using inflated radii) and orbit eccentricity (see Leinhardt 
& Richardson 2005). As was shown by Leinhardt et al. 
(2015), there is never a significant amount of mass in this 
unresolved component (see also section 3.2.2). 

In this work we track both the total mass of debris in 
each annulus, and the mass of core material (since this 
proved to be a problem in Paper I). Additionally, each 
particle and debris annulus has an ‘origin histogram’, 
that tracks the fraction of that object’s mass originat¬ 
ing from within each 0.1 AU annulus (see Leinhardt & 
Richardson 2005 for details). 

2.4. Mantle stripping law 

Differentiated planetesimals are modelled with core 
and mantle components, which are tracked indepen¬ 
dently. In a collision between two such planetesimals core 
and mantle material may be exchanged, or distributed 
amongst fragments. Marcus et al. (2010) analysed SPH 
collisions between differentiated bodies that started with 
separate core and mantle components, tracking the two 
types of particles to determine the result of the collision. 
They found that the fractional mass of the iron core of 
the largest post-collision remnant scaled with impact en¬ 
ergy (that is, the more energetic the impact, the higher 
the core fraction of the largest remnant), and lies be¬ 
tween two limiting cases: model 1, in which cores always 
merge; and model 2, in which core material from the 
projectile is only accreted if the largest remnant is more 
massive than the original target (see Marcus et al. 2010; 
Paper I). As in Paper I we assume that any core ma¬ 
terial not accreted by the largest remnant is distributed 
to the smaller collision remnants (ensuring mass is con¬ 
served), starting by placing as much as possible in the 
second largest. 

In Paper I the two empirical models from Marcus et al. 
(2010) describing the mass of the iron core of the largest 
remnant of a collision were applied to each collision in the 
simulation in a post-processing step. This proved to be 
problematic as it could not account for the composition 
of the unresolved material^; this had to be estimated in 
order to calculate embryo compositions. To avoid this we 
have incorporated the models used in Paper I into our 

^ The simulations conducted for Paper I had no compositional 
information at runtime, and so there was no composition for the 
unresolved material which was accreted continuously between col¬ 
lisions. The post-processing could determine the compositional 
change due to collisions, but the fraction of mantle or core mate¬ 
rial in the accreted unresolved material was unknowable. 


A/'-body code, which can now track the core mass of both 
resolved particles and debris annuli directly, applying the 
mantle stripping model as each collision occurs. 

In Paper I very little difference was found between the 
final core fractions of embryos using the two models from 
Marcus et al. (2010). We investigate the difference again 
in this work, and also define a third model as the mean 
of the original two, Mcore,lr ,3 ~ (^core,lr,l A[core,lr,2)/2 
(see Table 3). 

Collisions are not the only process that can effect bulk 
composition and core-mantle ratio of planetesimals; the 
initial oxidation state may have varied across the terres¬ 
trial planet feeding zone, but we do not know what any 
such distribution would have been. For this reason, and 
so that we will see the effects due to collisions alone, we 
assume the same initial core fraction across the terrestrial 
planet region. We explore two cases for this initial core 
fraction (both for planetesimals and unresolved debris) 
matching the values used in Paper I: 0.35, representative 
of high iron content chondrites (EH), and 0.22, represen¬ 
tative of low iron content chondrites (LL). 

Core and mantle material are always isolated in our 
model, and we do not allow any mixing or reequilibration 
to occur. Accreted core material is assumed to merge in¬ 
stantly with the core of the target since we expect planets 
to be (at least partially) molten during accretion. This 
is a significant simplification for some proposed scenar¬ 
ios (in which ongoing equilibration between metal and 
silicate alters the metal-silicate ratios, e.g. Rubie et al. 
2015), but a reasonable one given that the timescales in¬ 
volved and the degree to which mixing would occur are 
still poorly understood (e.g. Dahl & Stevenson 2010). 

2.5. The Grand Taek 

Walsh et al. (2011) presented the Grand Tack model of 
the evolution of the Solar System in which the inner plan- 
etesimal disc is truncated by the inwards-then-outwards 
migration of Jupiter and Saturn. In this scenario, Jupiter 
forms first, creating a gap in the disc allowing it to mi¬ 
grate inwards via type II migration. Saturn grew more 
slowly and later began migrating inwards more rapidly 
than Jupiter, allowing Saturn to ‘catch up’ with Jupiter, 
becoming trapped in the 2:3 mean motion resonance. 
This occurs as Jupiter reaches 1.5 AU; the gaps in the 
gas disc opened by the gas giants then overlap, chang¬ 
ing the balance of the torques, and causing the migra¬ 
tion direction to reverse. Jupiter and Saturn then mi¬ 
grate outwards together as the gas disc dissipates, leav¬ 
ing the gas giants close to their present day orbits, and 
the outer Solar System in the appropriate configuration 
for a late instability (associated with the late heavy bom¬ 
bardment). Jupiter migrating inwards to 1.5 AU before 
‘tacking’ truncates the inner disc at ^lAU, producing 
the much sought after small Mars (Walsh et al. 2011; 
O’Brien et al. 2014). 

In the simulations that employ the Grand Tack we in¬ 
clude Jupiter as an Wbody particle and model Jupiter’s 
migration by giving it an additional acceleration to 
match the velocity changes described by Walsh et al. 
(2011). Unlike previous Grand Tack simulations e.g. 
Walsh et al. (2011), O’Brien et al. (2014), Rubie et al. 
(2015), which begin with a possibly unrealistic bi- 
modal distribution of planetesimals and embryos and al¬ 
low Jupiter to begin migrating at the start, we begin 
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Table 3 

Mantle stripping models from Marcus et al. (2010). Model 3 is defined as the average of the original two: 

^core,lr,3 ~ core,lr,l + Mcore,lr, 2 )/ 2 . 


Collision outcome 

Mcorejr Model 1 

Mcorejr Model 2 

Perfect merging 

Partial accretion 

Hit & run 

Erosive, supercatastrophic disruption 

Mcore,targ + Mcore,proj 
min(Mij., Mcore,targ + Mcore,proj) 
Mcore,targ 

min(Mij., Mcore,targ + Mcore,pro.i) 

Mcore,targ + Mcore,proj 

Mcore,targ + min(Mcore,proj 5 - Mcore,targ) 

Mcore,targ 

min(Mij .5 Mcore,targ) 


with a less evolved planetesimal disc, and allow embryos 
to begin growing via collisions while Jupiter remains at 
3.5 AU for several Myr before it begins migrating. 

In this work we test the effect of Jupiter’s inward then 
outward migration starting at two different times - since 
Jupiter’s formation time, and the gas disc lifetime are 
uncertain - the commonly quoted time of 3 Myr (56k yr 
simulation time), and a date’ Grand Tack in which the 
migration begins 9 Myr (220 kyr simulation time) after 
the birth of the solar nebula. We also test two timescales 
for migration, both the standard 100 kyr (2.8 kyr simu¬ 
lation time) value of Walsh et al. (2011), and a ‘slow’ 
migration timescale of 600kyr (17kyr simulation time). 
Using the expansion factor method makes matching the 
timescales problematic, we have opted to accelerate the 
migration by the same factor by which the accretion 
timescale is reduced. Comparison of the slow and stan¬ 
dard timescales will show the effect of slower migration, 
but should also reveal any effect of the accelerated evo¬ 
lution relative to the migration caused by the expansion 
factor (the slow timescale is equivalent to the nominal 
migration rate accelerated only by a factor /). 

The acceleration, Uin(out): given to Jupiter to model its 
inward (outward) migration is described by: 

'^3.5AU — '^1.5AU /o^ 


T 

where ui.sau is Jupiter’s Keplerian velocity at 1.5 AU 
- the location at which it ‘tacks’, t is the time since 
the start of the simulation, ttack is the time at which 
Jupiter’s migration reverses direction, and r is the migra¬ 
tion timescale. Jupiter begins the simulation at 3.5 AU, 
reaches an innermost semi-major axis of 1.5 AU, and ends 
the simulation after the Grand Tack at ^5.2 AU. Jupiter 
feels the gravitational influence of the planetesimals just 
as any planetesimal or embryo, however, any object that 
collides with Jupiter is forced to merge. 

We do not include Saturn nor the ice giants in this 
work as we expect them to have no significant effect on 
the terrestrial planet region, as discussed by Walsh et al. 
(2011). We also ignore the outer planetesimals originally 
situated beyond the orbit of Jupiter, which would be 
numerically very expensive to model at these resolutions, 
as it is not expected that a significant number would 
be scattered into the terrestrial region (^3 percent of 
final planet mass, O’Brien et al. 2014; Ruble et al. 2015), 
and most would be accreted after the end time of our 
simulations. 


2.6. Calm disc 

In the calm disc scenario the initial planetesimal disc 
extends from 0.5-1.5 AU, with no giant planets included. 


2.7. Gas drag 

In this paper we present some simulations that include 
a gas density distribution within the protoplanetary disc 
(the effects of gas drag were not included in Paper I). 
In these simulations we calculate the drag force on the 
planetesimals following the prescription of Adachi et al. 
(1976) and Brasser et al. (2007). 

The drag force is applied to each body in the simulation 
as an additional acceleration in the direction opposite to 
that of the body’s relative motion through the gas, 

*D = - Vgf (4) 

where m and r are the planetesimal mass and radius, Cr> 
is the dimensionless drag coefficient, pg is the gas density, 
and V — Vg is the velocity of the planetesimal relative to 
the gas velocity. The drag acceleration is calculated using 
the un-inflated planetesimal radius, and is applied only 
to resolved particles; unresolved material is assumed to 
accrete before its orbit is significantly affected by drag, 
or to be replenished by inward drift from larger radii. 

The gas profile we adopt has one of two forms, either 
a power law matching the surface density of the solar 
nebula, or a third order polynomial obtained from a fit 
to the profile used by Walsh et al. (2011), that matches 
the results of hydrodynamic simulations from Morbidelli 
& Grida (2007). The MMSN gas density is given by 
(Brasser et al. 2007), 


pg = pg,is ^ 


( 5 ) 


where pgp is the initial gas density at 1 AU, s is the 
distance from the central star in the (x, ^)-plane, and 
2^5 = (0.047 5^/^) is the scale height of the gas disc, and 
the hydro dynamically derived gas density is given by. 




2 / 2 

s \ 


+ 25420— + 663.3 , ^ , 

5J J y/'KZs 


( 6 ) 


if s <= 0.9473 5j, or. 


5.2 e ^ 

Pg = — 90 — 

<sj \/'KZs 


( 7 ) 


otherwise, where sj is the instantaneous s of Jupiter. 
These equations provide a good approximation for the 
gas profile given in Walsh et al. (2011) interior to the 
orbit of Jupiter, the region with which we are concerned. 
We adopt a temperature profile T = 280 5“^/^ K as in 
Brasser et al. (2007), and thus the pressure gradient is 
given by 7^ = 2.13 X 10“^ The gas velocity is then 
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obtained from, 

Vg = VK \/l - 27], (8) 

where Vk is the Keplerian velocity at a distance s from 
the star. 

In our Grand Tack simulations the gas density be¬ 
gins to decay exponentially as Jupiter migrates outwards 
through the disc, whereas in the calm disc simulations 
the gas density either decays over the same time period 
or remains constant throughout. 

3. RESULTS 

The simulations discussed in the rest of this paper are 
summarised in Table 4. The larger disc and added com¬ 
plexity required for the Grand Tack simulations leads to 
substantially longer running times, and so we were only 
able to conduct a small number of high resolution ver¬ 
sions of these simulations. The similarity between the 
results of high and low resolution runs both for the calm 
disc and the Grand Tack model allows us to be confi¬ 
dent that the low resolution simulations are sufficient to 
demonstrate the behaviour of the system. 

3.1. Planetary embryo evolution 

Figure 2 shows the semi-major axis vs eccentricity at 
five snapshots in time of a system of planetesimals evolv¬ 
ing in the Grand Tack scenario. This broadly reflects 
the results of the Grand Tack as seen in other works, 
e.g. Walsh et al. (2011), O’Brien et al. (2014), namely, 
Jupiter truncates the inner disc and scatters some ma¬ 
terial outwards, leaving some higher eccentricity plan¬ 
etesimals between 2-3 AU to populate the asteroid belt. 
Note that the banded structure seen beyond ^1.3 AU in 
the third panel of Figure 2, and the grouping of excited 
planetesimals seen at ^1.9 AU in the bottom two panels 
are caused by orbital resonances with Jupiter during its 
(outward) migration, and would likely be disrupted by 
the presence of Saturn. 

Gomparing the evolution under the Grand Tack model 
(Figure 2) to the evolution of a calm disc environment 
(Figure 3), it can be seen that planetesimal eccentric¬ 
ities become much higher in the Grand Tack scenario, 
and larger embryos grow more rapidly. Examining the 
colours of the planetesimals and embryos, which indi¬ 
cate the initial location of material, it is clear that some 
blue or green material originating beyond 1.5 AU has 
been pushed into the inner disc by Jupiter’s migration, 
and that the disc is more mixed compared to the calm 
disc scenario, with only orange, yellow or green objects 
present in the inner disc at the end of the Grand Tack 
simulation, compared to the calm disc in which the final 
colours of objects approximately match the initial colours 
at those locations. It is worth noting that this mixing 
may be reduced compared to what would be seen with 
the nominal migration rate as the efficiency of resonant 
transport may have been changed by the increase in 
the migration rate. 

The results of our simulations for calm discs (simu¬ 
lations 1-19) are broadly similar to previous investiga¬ 
tions e.g. Kokubo & Ida (2002); Leinhardt & Richard¬ 
son (2005), with the same general differences noted in 
Paper I, namely that the larger number of particles 
(smaller starting sizes) and the inclusion of a sophisti¬ 
cated collisional model with fragmentation increase the 


timescales for runaway and oligarchic growth, as would 
be expected. The main obvious difference, before we con¬ 
sider the planetesimal compositions, between the simu¬ 
lations presented here and those presented in Paper I is 
the larger number of resolved fragments. This is due to 
the improved implementation of fragmentation described 
by Leinhardt et al. (2015). 

The mass distributions at the end of simulations 21 
and 5 are shown in Figures 4 and 1. Again this shows 
that some embryos grow larger in the same time period 
under the Grand Tack scenario. Here we are seeing the 
effect of shepherding significant amounts of mass from 
beyond 1.5 AU into the inner region of the disc; that 
the Grand Tack simulations result in embryos that are 
already one Earth-mass or more in only ^20 Myr of evo¬ 
lution is the result of this higher mass in the inner disc. 
The distinct lines of dark blue coloured planetesimals 
seen in the left hand panel of Eigure 4 show objects that 
underwent very little evolution in the outer part of the 
modelled disc (beyond ^2.5 AU) before being scattered 
onto highly eccentric orbits by Jupiter. They have subse¬ 
quently undergone very little or no interaction with other 
planetesimals, but could be responsible for late delivery 
of volatile elements (e.g. H, G). 

The right hand panel of Eigure 4 shows a very similar 
result to previous work investigating calm protoplanetary 
discs (e.g. Leinhardt et al. 2015; Paper I); embryos grow 
and detach from the mass distribution of the planetesi¬ 
mal population. As noted by Leinhardt et al. (2015) the 
embryos clean out the planetesimal population, starting 
in the inner disc where the evolution is faster, leading to 
an “embryo front” that propagates outwards through the 
disc, and there is no period where the total mass is split 
uniformly between planetesimals and embryos across the 
entire disc. This “inside-out” growth that results from 
the more primitive and realistic initial condition is in 
contrast to many recent A'-body models (e.g. Ghambers 
2013; Jacobson & Morbidelli 2014). This is also true in 
the Grand Tack scenario, but with the added complica¬ 
tion of the influx of material shepherded by Jupiter. 

3.2. Core fraetions 

Eigure 5 shows the fraction of each type of collision 
that occurred as a function of simulation time for the 
same Grand Tack simulation, 21, and calm disc simula¬ 
tion, 5 (both high resolution simulations that ignore gas 
drag, the effects of gas drag will be discussed below). It 
is immediately apparent that there are more erosive col¬ 
lisions (yellow) and supercatastrophic disruptions (red) 
under the Grand Tack scenario compared to the calm 
disc case. The left hand panel of Eigure 5 shows a spike 
in the number of collisions and fractions of supercatas¬ 
trophic disruptions and erosive collisions during Jupiter’s 
inward migration, and a higher fraction of these disrup¬ 
tive collision types than in the calm disc throughout the 
rest of the simulation. 

The collision history for the calm disc simulation shown 
in the right hand panel of Eigure 5 is very similar to that 
shown in Paper 1. As was seen in both Paper I and Lein¬ 
hardt et al. (2015) the fraction of disruptive collisions 
(yellow and red) increases with time as embryos grow, 
increasing the scattering of planetesimals. Note that the 
number of collisions decreases strongly with time in both 
scenarios (white lines), as the number of resolved parti- 
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a(AU) 

Figure 2. Evolution of the eccentricity vs semi-major axis for growing terrestrial planet embryos in the Grand Tack model (simulation 
21). The colours of the planetesimals represent the initial location of the material from which they are formed, as indicated by the colour 
bar. The sizes are proportional to the radii of the planetesimals, and the black centres represent their core mass fractions. The large 
black circle represents Jupiter, which is not to scale with the rest of the objects. Black and red labels show the simulation and effective 
timescales. The inset in the first panel shows a close-up of the inner part of the disc. (An animation of this figure is available online.) 
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Table 4 

Summary of simulations presented in this work. 


Simulation 

Ni 

Initial core fraction 

Mantle stripping 

Gas disc 

Jupiter 

1 

10 000 

0.22 

1 

None 

- 

2 

10 000 

0.22 

2 

None 

- 

3 

10 000 

0.22 

3 

None 

- 

4 

10 000 

0.35 

1 

None 

- 

5 

100 000 

0.22 

3 

None 

- 

6 

100 000 

0.22 

3 

None 

- 

7 

100 000 

0.22 

3 

None 

- 

8 

100 000 

0.35 

3 

None 

- 

9 

100 000 

0.35 

1 

None 

- 

10 

100 000 

0.35 

2 

None 

- 

11 

100 000 

0.35 

3 

None 

- 

12 

10 000 

0.22 

3 

Constant MMSN 

- 

13 

10 000 

0.35 

3 

Constant MMSN 

- 

14 

100 000 

0.22 

3 

Constant MMSN 

- 

15 

100 000 

0.35 

3 

Constant MMSN 

- 

16 

10 000 

0.22 

3 

Dissipating MMSN 

- 

17 

10 000 

0.35 

3 

Dissipating MMSN 

- 

18 

100 000 

0.22 

3 

Dissipating MMSN 

- 

19 

100 000 

0.35 

3 

Dissipating MMSN 

- 

20 

10 000 

0.22 

3 

None 

Fast, early GT 

21 

100 000 

0.22 

3 

None 

Fast, early GT 

22 

10 000 

0.22 

3 

None 

Slow, early GT 

23 

100 000 

0.22 

3 

None 

Slow, early GT 

24 

10 000 

0.22 

3 

None 

Fast, late GT 

25 

10 000 

0.22 

3 

MMSN 

Fast, early GT 

26 

10 000 

0.22 

3 

Hydrodynamic 

Fast, early GT 

27 

10 000 

0.22 

3 

Hydrodynamic 

Fast, early GT 

28 

10 000 

0.35 

3 

Hydrodynamic 

Fast, early GT 

29 

10 000 

0.22 

3 

Hydrodynamic 

Slow, early GT 

30 

10 000 

0.22 

3 

Hydrodynamic 

Fast, late GT 
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Figure 3. Evolution of the eccentricity vs semi-major axis for growing terrestrial planet embryos in a calm protoplanetary disc (simulation 
5). The colours of the planetesimals represent the initial location of the material from which they are formed, as indicated by the colour 
bar. The sizes are proportional to the radii of the planetesimals, and the black centres represent their core mass fractions. Black and red 
labels show the simulation and effective timescales. The inset in the first panel shows a close-up of the inner part of the disc. Note that 
the range of semi-major axis shown here is smaller than that in Figure 2, but the colour scales for this region of the disc are the same. (An 
animation of this figure is available online.) 
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Figure 4. Mass distribution of planetesimals and embryos against semi-major axis at the end (600 000 yr simulation time) of a Grand 
Tack simulation (21, left) and a calm disc simulation (5, right). The sizes of points are proportional to the radii of the planetesimals, and 
their colour indicates the mass-weighted radial source composition. (Animations of these figures are available online.) 
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Figure 5. Collision history for a fast, early Grand Tack simulation (21, left) and a calm disc simulation (5, right) as a function of simulation 
time. The coloured histogram shows the percentage (left axes) of each type of collision occurring in each 20kyr time bin (simulation time). 
The white histogram shows the number of collisions that occurred in each time bin (right axes). The white dashed lines indicate the times 
of the start of Jupiter’s migration, the reversal of the migration direction, and the approximate end of the outward migration (after ten 
outward migration timescales). Note that the first two of these lines are very close together on this scale, and that both occur in the third 
histogram bin. 
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cles falls. Again the fraction of accretional collisions is 
approximately constant with time. 

Unsurprisingly the dynamical excitation of the Grand 
Tack leads to a much more violent environment. The sig¬ 
nificant numbers of planetesimals with very large black 
cores visible in Figure 2 immediately suggests that migra¬ 
tion of Jupiter has had a major effect on the planetesimal 
core fractions. The influence of the migration on core 
fraction of the growing embryos is more clearly evident 
in Figure 6 . In this case with no gas and the standard 
fast, early migration (simulation 21 , left panel of Figure 
6 ), there is a large variation in core fraction following 
Jupiter’s migration (blue triangles), and some evolution 
towards lower core fractions is noticeable as the simu¬ 
lation proceeds (black circles). Even though the largest 
core fraction of an embryo falls during this later evolu¬ 
tion, the maximum is still higher than is ever reached in 
the calm disc, and the range in embryo core fraction is 
substantially larger (0.15-0.25 compared to 0.20-0.23 for 
the calm disc). 

Figure 7 shows the final core fractions of all resolved 
objects against their final semi-major axis, with the size 
of the points representing the radii of the planetesi- 
mals/embryos. In both simulations there is a small scat¬ 
ter in final core fractions of the embryos around the ini¬ 
tial value of 0.22 (as seen in Figure 6 ), with the largest 
bodies in the left hand panel having core fractions very 
close to the initial value. The smallest resolved planetes¬ 
imals, however, show a huge variation in core fraction, 
with some planetesimals composed entirely of iron, and 
some entirely of silicate mantle material (again in the 
left hand panel there are many largely unevolved blue 
planetesimals with small semi-major axis on high eccen¬ 
tricity orbits that still have approximately their initial 
core fraction). 

3.2.1. The mantle stripping models 

We found, as in Paper I, that the choice of mantle 
stripping law made little difference to the final core frac¬ 
tions of embryos. Since the original two models from 
Marcus et al. (2010) should over and underestimate the 
core fraction of the largest remnant, all later simulations, 
and those discussed below, use the average model, model 
3 (see Table 4). 

3.2.2. Unresolved debris 

Figure 8 shows the mass of unresolved debris present in 
the inner disc (0.5-1.5AU) as the simulations progress. 
There is a clear spike in the left panel of Figure 8 corre¬ 
sponding to Jupiter’s migration, and a number of smaller 
spikes in both simulations caused by collisions. 

In both simulations shown in Figure 8 the initial core 
fraction of all planetesimals and the initial unresolved 
material is 0 . 22 ; the average core mass fraction of the 
debris after 100 kyr is 0.14 and 0.10 in the Grand Tack 
and calm disc simulations respectively, and the average 
mass in unresolved debris after 100 kyr is 4x 10 “^ M 0 and 
2x10“^ M 0 . The unresolved material is enriched in sili¬ 
cate mantle material, which is as we would expect since 
it is mantle that is lost first in fragmenting collisions. As 
was shown by (Leinhardt et al. 2015), the total mass in 
the unresolved component is always a small fraction of 
the total system mass, and never builds to a significant 
mass. 


3.3. Migration timescales 

Figures 9 and 10 show the collision histories and em¬ 
bryo core fractions resulting from slower migration (sim¬ 
ulation 23) and late migration (simulation 24). 

The left hand panel of Figure 9 shows that the slower 
migration leads to a larger spike in the fraction of disrup¬ 
tive collisions (yellow and red), and a decrease in the frac¬ 
tion of ‘perfect’ hit and run collisions (light blue) com¬ 
pared to the standard Grand Tack scenario (left panel 
of Figure 5). The hit and run collisions in which the 
projectile remains intact (light blue) cannot affect the 
composition of either body in our model, whereas the 
partial accretion collisions (dark blue) that replace them 
in this slower migration scenario can. The slower migra¬ 
tion leads to a noticeable increase in the maximum final 
core fraction (0.26-0.28 for slow migration compared to 
0.24-0.25 for the fast scenario), as shown in the left panel 
of Figure 10 compared to the left panel of Figure 6 . 

The late migration (shown in the right panels of Fig¬ 
ures 9 and 10) causes an increase in the fraction of dis¬ 
ruptive collisions similar to the standard fast, early mi¬ 
gration simulations, but the embryo core fractions in this 
much more evolved planetesimal system do not undergo 
the same changes. The embryos have grown large enough 
during their longer period of evolution that the migra¬ 
tion has little effect on their composition. The variation 
in the final core fractions is significantly smaller than in 
the early Grand Tack simulations, and there is no large 
variation after Jupiter’s migration (green diamonds) to 
be averaged out by the end of the simulation (black cir¬ 
cles). 

3.4. The effects of gas drag 
3.4.1. Calm dise 

We have tested the influence of gas drag from a MMSN 
gas disc on the evolution of a calm planetesimal disc un¬ 
der two scenarios. In the first, the gas has a density that 
is constant in time throughout the duration of the sim¬ 
ulations; in the second the gas begins to dissipate after 
2.1 Myr (58 kyr simulation time) - approximately 3Myr 
after the birth of the solar nebula - the surface den¬ 
sity then decays exponentially on a timescale of 100 kyr 
(2.8kyr simulation time). These timescales were chosen 
to match those we use in the Grand Tack simulations 
(see section 3.4.2), which match those from Walsh et al. 
( 2011 ). 

The gas lifetime in the first scenario is probably unre¬ 
alistically long compared to observed gas disc lifetimes 
(^ 6 Myr, Haisch et al. 2001), but provides a useful in¬ 
sight into the effects of gas drag on disc and planetary 
embryo evolution. 

This scenario produces smaller oligarchs compared to 
the case without gas. Such an effect of gas drag was 
discussed by Tanaka & Ida (1997), it slows the growth 
of embryos as planetesimals become trapped outside the 
embryo feeding zones. The left hand panel of Figure 11 
shows that the inner region of the disc has been cleared 
to a greater degree than in simulations without gas drag 
(right panel of Figure 7). This reflects the more rapid 
evolution of the system due to gas drag, as seen by 
Wetherill & Stewart (1993) and Kokubo & Ida (2000). 

Embryos grow more quickly in the runaway growth 
phase due to the eccentricity damping caused by the gas. 
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Figure 6. Core fractions of embryos at four times during a Grand Tack simulation (21, left) and a calm disc simulation (5, right). Note 
that the red squares in the left hand panel represent a time before the migration begins. The grey area shows the final total range for the 
embryos in this region - defined as objects with a mass of at least 2 lunar masses. In the simulations depicted in this figure the initial core 
fraction was set to 0.22, indicated by the solid line. 



a(AU) 


Figure 7. Core fractions of planetesimals and embryos at the end (600 000 yr simulation time) of a Grand Tack simulation (21, left) 
and a calm disc simulation (5, right). The sizes of points are proportional to the radii of the planetesimals, and their colour indicates 
the mass-weighted radial source composition. The solid line indicates the initial core fraction. (Animations of these figures are available 
online.) 


however, during the oligarchic growth phase the embryos 
do not grow to be as massive as they are unable to sweep 
up the planetesimals. We are thus left at the end of our 
simulations with a system that has features that are in¬ 
dicative of both more and less evolved systems than the 
case without gas. 

Figure 12 shows the significant change in the fraction of 
collision types when gas is included; the left hand panel 
shows that gas drag results in significantly fewer erosive 
and disruptive collisions (yellow and red), and a greater 
fraction of perfect merging collisions (black) compared 
to the case without gas drag (Figure 5, right). Gas drag 
damps the eccentricities and inclinations of the planetes¬ 
imals, resulting in lower average collision velocities. The 
lack of these erosive and disruptive collisions in the con¬ 
stant gas scenario results in very few collisions that can 
significantly affect the core-mantle balance of the grow¬ 


ing embryos. The final distribution of core fractions (left 
panel of Figure 11) therefore shows very little change 
from the initial value, and much less variation than the 
simulations without gas drag (Figure 7, right panel). 

The second scenario, in which the gas dissipates early 
in the simulation, unsurprisingly has similarities to both 
the constant gas scenario and the ‘no gas drag’ scenario. 
The gas drag again facilitates the rapid growth of small 
embryos, however, once the gas density begins to decay 
the embryo eccentricities can be excited, allowing them 
to collide, producing a smaller number of larger embryos 
that go on to continue clearing the planetesimal popula¬ 
tion. 

The collision histogram (right panel of Figure 12) looks 
similar to the no gas drag case (Figure 5, right panel) 
once the gas begins to dissipate, but in the first 58kyr 
(2.1 Myr effective time) - during which the majority of 
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Figure 8. Unresolved debris mass as a function of simulation time for a Grand Tack simulation (21, left) and a calm disc simulation (5, 
right). The black line indicates the total mass in unresolved debris between 0.5 and 1.5 AU, and the red line indicates the corresponding 
mass of core material. 
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Figure 9. Collision history of Grand Tack simulations as a function of simulation time, for a high resolution simulation with slow, early 
migration (23, left), and a low resolution simulation with fast, late migration (24, right). The coloured histogram shows the percentage 
(left axes) of each type of collision occurring in each 20kyr time bin (simulation time). The white histogram shows the number of collisions 
that occurred in each time bin (right axes). The white dashed lines indicate the times of the start of Jupiter’s migration, the reversal of the 
migration direction, and the approximate end of the outward migration (after ten outward migration timescales). Note that low resolution 
simulations tend to have more perfect hit and run collisions due to the smaller range of particle sizes. 


the collisions occur - there are fewer erosive events and 
more perfect merging collisions. This likely leads to the 
smaller variation in final core fraction shown in Figure 
11 (right hand panel) for this dissipating gas scenario 
compared to the case without gas drag (Figure 7, right 
panel). 

3.4.2. The Grand Tack scenario 

Including gas in the standard fast, early migration sce¬ 
nario leads to embryo core fractions that evolve in a 
similar way to the simulations with slower migration, 
with the larger variation and enhanced final core frac¬ 
tions compared to the case without gas drag (see Figure 
13). Whilst the collision histograms look most similar 
to the fast, early migration without gas, the inclusion of 
gas drag has prevented the enhanced core fractions from 
being reduced by the end of the simulation. The results 
are similar for both a MMSN gas density profile (sim¬ 


ulation 25) and the hydrodynamic derived gas density 
profile (simulation 26) used by Walsh et al. (2011). 

Including gas drag in the Grand Tack scenario in¬ 
creases the maximum core fraction compared to the sce¬ 
nario without gas drag; this is the opposite effect to that 
in the calm disc scenario, in which the damping of eccen¬ 
tricities caused by the drag force leads to a smaller varia¬ 
tion in core fraction. The excitation caused by Jupiter’s 
migration has a substantially larger effect than the damp¬ 
ing caused by the drag force in a MMSN gas disc. 

4. DISGUSSION 

The mean and range of embryo core mass fractions 
from all of the simulations are shown in Figure 13 (where 
we define embryos as bodies with a mass of at least 2 
lunar masses; note that this mean value is not weighted 
by the embryo mass, and that both core material and 
mantle material are conserved in these simulations). It is 
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Figure 10. Evolution of the core fractions of embryos during Grand Tack simulations for a high resolution simulation with slow, early 
migration (23, left), and a low resolution simulation with fast, late migration (24, right). The grey area shows the final total range for the 
embryos in this region - defined as objects with a mass of at least 2 lunar masses. In the simulations depicted in this figure the initial core 
fraction was set to 0.22, indicated by the solid line. (Animations of these figures are available online.) 
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Figure 11. Core fractions of planetesimals and embryos at the end (600 000 yr simulation time) of calm disc simulations with constant gas 
(14, left), and dissipating gas (18, right). The sizes of points are proportional to the radii of the planetesimals, and their colour indicates 
the mass-weighted radial source composition. The solid line indicates the initial core fraction. (Animations of these figures are available 
online.) 


important to note that in all cases we find both embryos 
that show an enhanced core fraction, and embryos that 
show a decreased core fraction. 

It is clear from Figure 13 that the dynamical excitation 
caused by Jupiter’s Grand Tack leads to a higher maxi¬ 
mum and larger range of embryo core fractions compared 
to the calm scenario. This is not entirely surprising, 
Jupiter’s Grand Tack excites the orbits of many plan¬ 
etesimals and embryos, which makes them more likely 
to undergo more energetic collisions. This leads to larger 
numbers of disruptive collisions, and it is these disruptive 
collisions that can strip off large amounts of mantle, and 
even excavate material from the cores of planetesimals 
or embryos, allowing core material to be redistributed to 
other objects. 

Dwyer et al. (2015) investigated the compositional evo¬ 
lution of terrestrial planets during the late stages of ac¬ 


cretion by post-processing low N simulations of the giant 
impact phase. In this work we have studied an earlier 
stage of planet formation: the runaway and oligarchic 
growth phases, so a direct comparison is not possible. 
Dwyer et al. determined the core fractions of bodies 
resulting from imperfect collisions using a similar ap¬ 
proach to our model 1, however, rather than using the 
unresolved debris method, the simulations they analysed 
either placed all the fragment mass into equal mass bod¬ 
ies, or, if this mass was below the resolution limit, it was 
simply accreted by the target (or projectile in the case 
of hit and run collisions, Ghambers 2013). This means 
that a significant amount of mantle had to be stripped 
before a change in composition would occur. They found 
a greater scatter in silicate mass fraction for lower mass 
bodies; similarly we find a greater variation in composi¬ 
tion for smaller bodies, though it should be noted that we 
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Figure 12. Collision history for high resolution calm disc simulations with constant gas (14, left), and dissipating gas (simulation 18, 
right), as a function of simulation time. The coloured histogram shows the percentage (left axes) of each type of collision occurring in each 
20kyr time bin (simulation time). The white histogram shows the number of collisions that occurred in each time bin (right axes). 
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Figure 13. Mean, minimum and maximum final embryo core 
mass fractions from all simulations. The dotted lines indicate the 
initial core fractions. Model 1, 2 and 3 refer to the mantle stripping 
laws. 
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cover an earlier stage of evolution here, and that our mass 
limit is more than 100 times smaller. In the Grand Tack 
model we find some embryos with masses approaching 
Earth-mass (although their evolution is not complete), 
which tend to show only a small variation in core frac¬ 
tion similar to the results of the calm disc simulations 
(e.g. Figure 7). The smaller embryos in the Grand Tack 
generally exhibit a larger range of core fractions, though 
these are still modest in comparison to the range shown 
by the smallest remaining planetesimals in both scenar¬ 
ios. 

The range of final embryo core fractions in the calm 
scenario is smaller than in Paper I, which was by ne¬ 


cessity an uncertain estimate due to the way in which 
accretion of unresolved material was handled. The up¬ 
dated range of final core fractions fit within the ranges 
from Paper I, but have more conservative maximum val¬ 
ues. Even with our direct tracking of the composition 
of the unresolved debris there are still several limitations 
with this method. The dynamics of this material and any 
collisional evolution are ignored. In a real system a small 
fraction of this mass would likely be ground down to dust 
and lost via Poynting-Robertson drag, and the composi¬ 
tion of individual sub-200-km bodies would also evolve. 
Since the unresolved material is enhanced in silicate com¬ 
pared to the initial conditions, any loss of this material 
from the system would further enhance core fractions of 
the resolved bodies. However, most of the mass in the 
unresolved debris would be in bodies with masses just be¬ 
low the resolution limit, only a tiny fraction of the mass 
would be in the form of dust that could be lost, and so 
we expect the effects of this to be minimal. The up¬ 
dated results for the calm scenario lead to smaller shifts 
in composition than seen in Paper I, but the collisional 
evolution of chondritic planetesimals could still produce 
embryos with Earth like composition in extreme cases 
with high initial core fractions. 

As has been noted, we see a much larger variation in 
core fraction in the Grand Tack models than in the calm 
scenario. The fact that the maximum is higher when 
Jupiter’s migration is slower, or when we include the ef¬ 
fects of gas drag, is more difficult to understand. Figure 
9 shows that there are fewer collisions after the Grand 
Tack in these cases compared to the scenario without 
gas drag, the slightly more evolved systems then undergo 
fewer collisions that could work to reduce the variation in 
core fraction that is present immediately after Jupiter’s 
migration in all four cases (Figures 6 and 10). In the 
case of slower migration, it is expected that the efficiency 
of resonant transport will be increased, which would lead 
to a larger number of excited planetesimals and hence a 
greater number of erosive collisions in the inner regions 
of the disc. Gas drag will cause the smallest fragments 
produced via collisions during Jupiter’s inward migration 
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to lose eccentricity and move inwards towards the Sun, 
and since these fragments will be dominated by man¬ 
tle material, gas drag can act to enhance the core frac¬ 
tions of embryos closest to the tack point. On the other 
hand, if the Grand Tack occurs late, too much evolution 
has already occurred, and the excitation is not sufficient 
to cause large variation in the core fractions of the now 
larger embryos. This suggests that the exact composition 
of the Earth could be sensitive to the timing of Jupiter’s 
migration under the Grand Tack scenario, and that with 
a very good understanding of the evolution of the Earth’s 
composition, it may be possible to ‘date’ the Grand Tack. 

Eigures 6 and 10 reveal a noticeable gradient in core 
fraction across the region occupied by the terrestrial 
planets, with lower core mass fractions interior to ^1 AU, 
and higher core fractions beyond this point. It is the 
region beyond ^lAU that is most affected by the mi¬ 
gration of Jupiter, with much of the material originally 
in this region excited onto high eccentricity orbits before 
relaxing as Jupiter migrates outwards again (see Eigure 
2). This causes some silicate rich fragments to be scat¬ 
tered out of this region (both inwards and outwards), 
leaving the disc enriched in iron between 1-1.5 AU, and 
silicate enriched interior to 1 AU. This is presumably sen¬ 
sitive to the location of the innermost point of Jupiter’s 
migration, which was specifically chosen by Walsh et al. 
(2011) to affect the disc outside 1 AU. It is also impor¬ 
tant to remember that we have used a uniform initial 
core fraction, and it is expected that the planetesimals 
in a real system would have some variation in initial core 
fraction (e.g. Righter et al. 2006). Our final solar system 
clearly shows that the effects of the giant impact stage 
are significant and stochastic, with the high Ee/Mg ratio 
for Mercury possibly related to a giant impact and the 
preserved embryonic Ee/Mg for Mars. 

4.1. Mixing and volatile delivery 

The ‘origin histogram’ (see Eigure 14) of each parti¬ 
cle in our simulations tracks the initial locations of the 
mass of which that particle is comprised, so we can easily 
examine our final embryos and determine from which re¬ 
gions of the disc they have accreted planetesimals. The 
first panel of Eigure 14 shows the composition in terms 
of radial location of embryos produced in calm disc sim¬ 
ulations. This shows a very similar result to that found 
by Leinhardt et al. (2015), that there is some localised 
mixing with most of the mass accreted by each embryo 
originating from the regions 0.1 AU either side of its po¬ 
sition at the end of the simulation. 

Eigure 14 (lower panel) shows the mixing in the same 
way for the Grand Tack model. It is clear that Jupiter’s 
migration in this model causes a much greater degree of 
mixing, with almost all embryos containing some signif¬ 
icant fraction of material originating beyond 2AU, and 
no pair of adjacent 0.1 AU regions dominating the mass 
of any embryo that remains in the terrestrial planet re¬ 
gion. Eigure 14 also suggests a much more uniform radial 
composition due to this increased mixing than is seen in 
any of the calm disc simulations. Again we note that the 
accelerated migration may have reduced the efficiency of 
transport and thus caused reduced mixing. 

In this work we have assumed a uniform initial core 
fraction across the planetesimal disc, this allows for a 
clearer investigation of the effects of collisional evolution. 


but in a real system we would expect some variation in 
oxidation state with heliocentric distance (e.g. Righter 
et al. 2006). This would be expected to have some effect 
on the final core fractions of the terrestrial planets, espe¬ 
cially in the calm scenario which shows little mixing. The 
much more mixed embryos produced in the Grand Tack 
simulations may not be affected significantly by this sim¬ 
plification, however, the effect on the resulting system of 
planets will depend on the nature of the initial variation. 

In the Grand Tack model all embryos acquire simi¬ 
lar fractions of ‘blue’ material (3-8 percent of final mass 
from beyond 2.5 AU, see Eigure 14), that is often as¬ 
sumed to contain more water than planetesimals formed 
closer to the Sun (e.g. Raymond et al. 2009). If plan¬ 
etesimals originating beyond 2.5 AU are 10 percent wa¬ 
ter by mass (as is commonly assumed in similar calcula¬ 
tions, e.g. O’Brien et al. 2006), the embryos at the end 
of these simulations would be 0.4-0.8 percent water, and 
the larger embryos would already have acquired enough 
to match the Earth’s water content. Here we have used 
only planetesimals originating interior to the orbit of 
Jupiter, O’Brien et al. (2014) similarly estimate that any 
Earth mass planet forming under the Grand Tack would 
acquire the required amount of water from planetesimals 
originating beyond the orbit of Jupiter (which we have 
neglected). However, it is uncertain what fraction of the 
volatiles from planetesimals accreted during this period 
would be retained by the final planet, and water from 
both populations of planetesimals may be required if a 
substantial fraction is lost due to collisions. 


4.2. The fate of planetesimals 

Bottke et al. (2006) proposed that the iron meteorites, 
which originate today from the inner asteroid belt, are 
remnants of planetesimals formed in the terrestrial planet 
region rather than fragments of bodies formed in the 
main belt. One outstanding question from Bottke et al.’s 
work is whether differentiated planetesimals from the ter¬ 
restrial planet region could actually produce these iron- 
rich bodies via disruption. The results of this work sug¬ 
gest that many iron-rich, and complimentary silicate- 
rich, fragments are indeed produced by the collisional 
evolution of the inner disc, as seen in Eigure 7. Though 
we see very few of these remnant planetesimals scattered 
into the asteroid belt region, it is important to note that 
the initial distribution of planetesimals in our calm disc 
simulations does not extend beyond 1.5 AU, and that 
such material may be important in scattering fragments 
into the main belt region (and keeping them there). 

It has recently been realised that many white dwarfs 
display extrasolar planetary debris in their spectra (e.g. 
Jura & Young 2014). This material shows a diversity 
in the iron contents of accreted planetesimal and small 
embryo sized bodies (e.g. Earihi et al. 2011; Gansicke 
et al. 2012; Raddi et al. 2015) that could represent left¬ 
over collisionally processed planetesimals similar to those 
produced in our simulations (e.g. Eigure 7). The pres¬ 
ence of differentiated rocky debris in white dwarf atmo¬ 
spheres appears to be a common phenomenon, which 
may suggest that many extrasolar planetary systems un¬ 
dergo similar collisional evolution to that shown in this 
work. 
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Semi-Major Axis (AU) 


Figure 14. Composition pies for each planetary embryo at the end point of a calm disc simulation without gas drag (5, upper), and a 
Grand Tack simulation with fast, early migration (21, lower). The embryos are shown at their instantaneous semi-major axis, with the 
size of each pie proportional to the object’s radius. Each sector of a pie shows the mass fraction of the embryo’s mantle or core originating 
from the region indicated by the colour bar. The histograms on the right show this breakdown for the indicated embryo, the dashed line 
showing its location. 


4.3. What does this mean for the bulk eomposition of 
planetary embryos? 

In order to estimate the bulk Fe/Mg ratios of the fi¬ 
nal planetary embryos, we use a series of starting bulk 
compositions for the planetesimals which span the range 
represented by chondritic meteorites. Following the tech¬ 
nique used in Paper I, for each bulk planetesimal compo¬ 
sition, a metallic core is formed of the proscribed mass 
(either 0.22 or 0.35, see section 2.4, with 10 percent of 
the core as light elements in the former or 5 percent light 
elements and 5 percent Si in the latter) and a residual 
mantle composition calculated. These compositions are 
used together with the modelled mass fractions of core 
and mantle in the final embryos (see Figure 13) to de¬ 
termine their bulk Fe/Mg. Finally, we apply the addi¬ 
tional constraint that the embryo must have enough Fe 
to match the Earth’s measured core mass fraction (the 
FeO/Fe ratio is adjusted in order to match the Earth’s 
core fraction of 32 percent, ensuring that the EeO content 
of the mantle is zero or greater). 

Eigure 15 illustrates the differences between starting, 
chondritic compositions (dashed lines) and bulk Earth 
composition (dark grey bar) indicating the need for some 
processing of chondritic material in order to produce the 
Earth. The ranges in Ee/Mg that result from a selection 
of our simulations are also shown in Eigure 15 and it is 
evident that collisional processing can produce embryos 
matching Earth’s Ee/Mg ratio. Earth-like Ee/Mg can 
be produced in extreme cases starting with planetesimal 
compositions close to bulk Earth in calm scenarios and 


more readily in the Grand Tack models. The collisional 
stripping of mantle caused by the Grand Tack can be 
sufficient to explain the non-chondritic bulk Ee/Mg ratio 
of the Earth starting with a wider range of planetesimal 
compositions. 

It is, however, important to note that the accretion 
of terrestrial planets is not complete at the end-point of 
these simulations, and that it is not necessarily the case 
that the embryos with enhanced core fractions represent 
the Earth. Gollisions during the giant impact phase may 
cause any compositional differences observed at the end 
of oligarchic growth to be averaged out again, or they 
may serve to enhance this variation (see Stewart & Lein- 
hardt 2012; Paper I for more detailed discussions of com¬ 
positional changes during the giant impact phase). 

The compositional variation that we see in these simu¬ 
lations also has strong implications for exoplanetary sys¬ 
tems. The bulk composition of planets could be signifi¬ 
cantly different from the expected composition based on 
the stellar photosphere. The GI chondrites are a good 
match for our Sun, but as Eigure 15 shows the bulk com¬ 
position of planetary embryos formed from Cl chondrites 
may not be the same. Given the interest in the compo¬ 
sitions of exoplanets and their atmospheres, it is impor¬ 
tant to consider the possible deviation from the compo¬ 
sition implied by their host stars, especially if they may 
have undergone dynamical processes similar to the Grand 
Tack. 
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Figure 15. Bulk Fe/Mg ratios for collisionally processed em¬ 
bryos. Using the range of final core fractions achieved in each 
simulation, the vertical bars show the range of possible Fe/Mg in 
embryos for a series of scenarios starting with different initial chon- 
dritic compositions (Cl through to EL), with the requirement that 
the embryo has sufficient Fe to match the Earth’s present core 
mass fraction. The dotted lines indicate the initial Fe/Mg for the 
five chondrite types, the solid black line represents bulk earth if 
the core contains 10 percent light elements, the dark grey shaded 
region shows the range for 0-15 percent light element content (as 
given by O’Neill & Palme 2008). If a range shown for a particu¬ 
lar simulation does not overlap with the dotted line indicating the 
initial Fe/Mg ratio for that chondrite it is not possible to reach 
the required 32 percent core fraction without stripping of man¬ 
tle. The upper limits always represent the Fe/Mg achieved from 
the embryo with maximum increase in core fraction, however the 
minimum values do not necessarily represent the embryo with the 
minimum core fraction from that simulation. The lower limit is 
set by the minimum final core fraction required, for that chon¬ 
drite type and initial oxidation, to retain sufficient Fe to match the 
Earth’s core fraction, with the requirement that this minimum is 
within the range of core fractions actually found. If there is no 
range shown for a simulation the maximum core fraction achieved 
was not sufficient to reach this minimum required Fe content. 

5. SUMMARY 

We have investigated the compositional evolution of 
growing terrestrial planet embryos through runaway and 
oligarchic growth both under the Grand Tack model 
(Walsh et al. 2011) and a model with no influence from 
giant planets. Our simulations conducted using PKD- 
GRAV (Richardson et al. 2000; Stadel 2001) with a state- 
of-the-art collision model (EDACM, Leinhardt & Stew¬ 
art 2012; Leinhardt et al. 2015) and mantle stripping 
laws (Marcus et al. 2010) show that planetary embryos, 
grown from differentiated planetesimals with an initially 
uniform core mass fraction, develop variations in their 
core fractions, and that both iron-rich and silicate-rich 
fragments are produced. 

In the Grand Tack scenario, and in extreme cases for 
the calm scenario, these variations in embryo core mass 
fraction can be sufficient to account for the Earth’s non- 
chondritic Ee/Mg ratio. 

In all of the simulations, we see both embryos with in¬ 
creased core fractions and embryos with decreased core 
fractions. Regardless of whether the compositional shift 
is in a direction that makes the Earth appear more chon- 
dritic, or less chondritic, it is clear that collisions could 
substantially alter the composition of planetary embryos 
from that of the initial differentiated planetesimals. 
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